function [theta, Sigma, Q_alpha, Q_sigma, logsigma, alpha] = SV_VAR_Simple(yraw,var_priors,Sigma,Q_alpha,Q_sigma, prior_SV,h,vardiag)
% SV_VAR produces a draw from the stochastic volatilities of a VAR
%   Detailed explanation goes here

%% Preliminaries

% Dimensions of Data and Model Objects

[Traw, n] = size(yraw);


if vardiag == 1 
    r = 0;
else
    r = (n*(n-1)/2);

end

y = yraw;

% Generate lagged Y matrix. This will be part of the X matrix
ylag = mlag2(y,2); % Y is [T x M]. ylag is [T x (Mp)]

X = ylag(2+1:Traw,:);

y = y(2+1:Traw,:); % This is the final Y matrix used for the VAR

T = Traw - 2;

A = nan(n,n,T);
Omega = nan(n,n,T);

alpha = nan(T,r);
logsigma = nan(T,n);

for t = 1:T
    
    Omegat = diag(diag(chol(Sigma(:,:,t),'lower')));
    At = Omegat/(chol(Sigma(:,:,t),'lower'));
        
    A(:,:,t) = At;
    Omega(:,:,t) = Omegat;

    if r>0 
        alpha(t,:) = At(tril(true(size(At)), -1))'; 
    end
    logsigma(t,:) = diag(log(Omegat))'; 
    
end

%% Priors 

    % Priors for starting values of TVP Parameters
    
    if r>0, alpha0 = zeros(r,1); end
    Sig_alpha0 = 10*eye(r);
    
    logsigma0 = -1;
    Sig_logsigma0 = eye(n);
    
    % Priors for Covariance Matrices
       
    sQ_alpha = prior_SV*eye(r);
    dfQ_alpha = 1;
    
    sQ_sigma = prior_SV*eye(n);
    dfQ_sigma = 1;


%% Start Gibbs Draw
    
    % Step 1 - Draw the Autoregressive Coefficients    
      
      factors = nan(T-h,n);
      for t = 1:T-h
      factors(t,:) = (inv(Omega(:,:,t))*A(:,:,t)*y(t,:)')';
      end
        
      [theta,~,~] = var_gibbs(factors,[],var_priors,2,eye(n),0);
      
      yhat = (y - X*theta);
      

    % Step 3 - Draw a history of off-diagonal elements, alpha
    
    if r>0
        % This part of the algorithm is based not on Primiceri 2005 but on
        % Section 3.3 of De Wind and Gambetti (2014)
        
        Z = nan(n,r,T);
        
            for t=1:T
                Astar = A(:,:,t) - eye(n);
                vecAstar = vec(Astar');
                Zt = -kron(eye(n),yhat(t,:));
                Z(:,:,t) =  Zt(:,abs(vecAstar)>1E-10);    % Transpose??           
            end
        
        [xtt, Sigtt,~,~,~] = KalmanFilter(yhat,[],[],alpha0,Sig_alpha0,[],eye(r),Z,[],Q_alpha,Omega,[]);
        alpha = Bai_Wang2(eye(r),Q_alpha,[],[],Sigtt,xtt,r);
        
    % Step 4 - Draw Q_alpha
    
    d_alpha = alpha-mlag2(alpha,1);
    d_alpha = d_alpha(5:end,:); % I supress a couple of observations to reduce dependency from initial conditions
    
    Q_alpha = DrawIW(d_alpha,sQ_alpha,dfQ_alpha);
    
    end
    
    % Step 5 - Draw the log-Standard Deviations
    
     ystar = nan(T,n);
     A = nan(n,n,T);

     for t=1:T
         At = eye(n);
         if r>0 
         At(tril(true(size(At)), -1)) = alpha(t,:) ;
         end
         A(:,:,t) = At;

         ystar(t,:) = yhat(t,:)*A(:,:,t)';
     end
     
         ystarstar = log(ystar.^2+1e-3);

    for i = 1:n    
    [S, logsigma(:,i)] = KimShephardChib(ystarstar(:,i),logsigma(:,i),logsigma0(i),Sig_logsigma0(i,i),Q_sigma(i,i));   
    end
%     logsigma = .5*logsigma;
         
    % Step 6 - Draw Q_sigma
    
    d_logsigma = logsigma-mlag2(logsigma,1);
    d_logsigma = d_logsigma(5:end,:);
    
    Q_sigma = DrawIW(d_logsigma,sQ_sigma,dfQ_sigma);
    
    % Step 7 - Reconstruct Sigma

        Sigma = nan(n,n,T);
        Omega = nan(n,n,T);

        for t=1:T
            Omega(:,:,t) = diag(exp(logsigma(t,:)));
            Sigma(:,:,t) = (inv(A(:,:,t))*Omega(:,:,t))*(inv(A(:,:,t))*Omega(:,:,t))';
        end
        
    % Saving and Plots



end

